
rm(list=ls())
setwd("~/Dropbox/Nigeria TfS Conjoint/Replication")

# list the packages that we load
# alphabetically for reproducibility
packages <- c('fastDummies', 'cjoint', 'monomvn', 'stringr', 
              'survey', 'cregg', 'estimatr', 'dplyr', 'coefplot', 'DT', 'MASS', 'haven','multiwayvcov','broom')
# call library on each package
purrr::walk(packages, library, character.only=TRUE)

# some packages we will reference without actually loading
# they are listed here for complete documentation
packagesColon <- c('knitr', 'magrittr', 'purrr', 'tibble', 'useful')

data<-read.csv("Tax_for_Services_Nigeria_final_analysis.csv")

data$whotaxpaid <- factor(data$whotaxpaid)
data$leveltax <- factor(data$leveltax)
data$deduct <- factor(data$deduct)
data$howtax <- factor(data$howtax)
data$healthcov <- factor(data$healthcov)
data <- within(data, healthcov <- relevel(healthcov, ref = 3))
data$perk <- factor(data$perk)

data<- subset(data, is.na(data$sampling_weight)==FALSE)

setwd("~/Dropbox/Nigeria TfS Conjoint/Results/")

data$rating <- as.numeric(gsub("\\.", "", data$rating))



####################################
############Figure 1 ###############
####################################

lm<-lm(rating ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data, weights=data$sampling_weight)
vcovCL<-cluster.vcov(lm, data$id)
coef<-coeftest(lm, vcovCL)

acme_estimates<-coef %>% tidy() %>% filter(term!="(Intercept)")%>% 
  mutate(lower=estimate-1.96*std.error) %>% mutate(upper=estimate+1.96*std.error) %>% 
  mutate(cate = substr(term, start = 1, stop = 3)) %>% mutate(feature_level=NA)

#  %>% filter(estimate!= "NA")
acme_estimates$feature_level[acme_estimates$term=="whotaxpaidYou would pay the tax to the Local Government Authority"]<-"Local Government Authority"                                                     
acme_estimates$feature_level[acme_estimates$term=="whotaxpaidYou would pay the tax to the State Internal Revenue Services"]<-"State Internal Revenue Services"                                              
acme_estimates$feature_level[acme_estimates$term=="leveltaxThe tax on income would be between N24000 and N36,000"]<-"N24000-N36000"                                           
acme_estimates$feature_level[acme_estimates$term=="leveltaxThe tax on income would be between N6000 and N12,000"]<-"N6000-N12000"                                                         
acme_estimates$feature_level[acme_estimates$term=="deductYou would receive low-cost healthcare, in which you co-pay N 750 each time you go to the doctor"]<-"Low-cost healthcare"                
acme_estimates$feature_level[acme_estimates$term=="howtaxYou would pay the tax to the bank"]<-"Bank"                                                                               
acme_estimates$feature_level[acme_estimates$term=="howtaxYou would pay via UUSD, which uses telecom service just like when you load credit on your cell phone"]<-"UUSD"            
acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare already 3 months before you start paying income tax"]<-"3 months before you pay"                             
#acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare at the moment you start paying income tax"]<-"at the moment start paying"                                       
acme_estimates$feature_level[acme_estimates$term=="perkWhen paying income tax, you would receive a state identity card that you could access new government services with"]<-"State identity card"
acme_estimates$feature_level[acme_estimates$term=="perkYour health coverage would be extended to one partner and four children"]<-"Health coverage extended to 1 partner and 4 children"
acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare 3 months after you start paying income tax"]<-"3 months after you pay"

acme_estimates <- acme_estimates %>% add_row(feature_level = "Who", estimate = NA, lower= 0, upper=0, cate="who")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Level", estimate = NA, lower= 0, upper=0, cate="lev")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Deduct", estimate = NA, lower= 0, upper=0, cate="ded")
acme_estimates <- acme_estimates %>% add_row(feature_level = "How", estimate = NA, lower= 0, upper=0, cate="how")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Health Coverage", estimate = NA, lower= 0, upper=0, cate="hea")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Perk", estimate = NA, lower= 0, upper=0, cate="per")

acme_estimates <- acme_estimates %>% add_row(feature_level = "Cooperative or association", estimate = 0, lower= 0, upper=0, cate="who")
acme_estimates <- acme_estimates %>% add_row(feature_level = "N12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Free healthcare", estimate = 0, lower= 0, upper=0, cate="ded")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Tax officer or collecting agent", estimate = 0, lower= 0, upper=0, cate="how")
acme_estimates <- acme_estimates %>% add_row(feature_level = "At the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Healthcare is kept for a certain time when income is lost", estimate = 0, lower= 0, upper=0, cate="per")


acme_estimates <- acme_estimates %>% mutate(feature_level = factor(feature_level,levels=c( 
  "State identity card", "Health coverage extended to 1 partner and 4 children", "Healthcare is kept for a certain time when income is lost","Perk",
  "3 months before you pay",  "3 months after you pay","At the moment start paying","Health Coverage",
  "Bank","UUSD", "Tax officer or collecting agent","How",
  "Low-cost healthcare","Free healthcare","Deduct",
  "N24000-N36000", "N6000-N12000" ,"N12000-N24000","Level",
  "Local Government Authority","State Internal Revenue Services","Cooperative or association","Who")))

acme_estimates %>% 
  #filter(estimate!=0) %>% 
  ggplot(aes(estimate, feature_level)) +
  geom_point(position=position_dodge(width = .4),size=4, alpha=0.8) +
  ggstance::geom_pointrangeh(aes(xmin=lower, xmax=upper), position=position_dodge(width = .4), size=0.7) +
  geom_vline(xintercept=0, linetype="dashed", color="grey") +
  theme(legend.title=element_blank(),
        legend.text=element_text(size=15),
        axis.title=element_text(size=15),
        axis.text=element_text(size=15))+
  xlab("")+
  ylab("")+
  scale_y_discrete(labels=rev(c(expression(paste(bold("Who:"))),
                                expression(paste("Cooperative or association")),
                                expression(paste("State Internal Revenue Services")),
                                expression(paste("Local Government Authority")),
                                expression(paste(bold("Level:"))),
                                expression(paste("N12000-N24000")),
                                expression(paste("N6000-N12000")),
                                expression(paste("N24000-N36000")),
                                expression(paste(bold("Deduct:"))),
                                expression(paste("Free healthcare")),
                                expression(paste("Low-cost healthcare")),
                                expression(paste(bold("How:"))),
                                expression(paste("Tax officer or collecting agent")),
                                expression(paste("UUSD")),
                                expression(paste("Bank")),
                                expression(paste(bold("Health Coverage:"))),
                                expression(paste("At the moment start paying")),
                                expression(paste("3 months after you pay")),
                                expression(paste( "3 months before you pay")),
                                expression(paste(bold("Perk:"))),
                                expression(paste("Healthcare is kept for a certain time when income is lost")),
                                expression(paste("Health coverage extended to 1 partner and 4 children")),
                                expression(paste("State identity card"))
  )))+
  guides(color=guide_legend(reverse = TRUE))-> p
p

####################################
############Table 2 ###############
####################################

acme_rounded <- acme_estimates

numeric_cols <- c("estimate", "std.error", "statistic", "p.value", "lower", "upper")
for(col in numeric_cols) {
  acme_rounded[[col]] <- round(acme_rounded[[col]], 2)
}

print(acme_rounded)


####################################
############Figure 2 ###############
####################################

acme_reg_old <- function(data_input){
  data<-data_input
  lm<-lm(rating ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data)
  vcovCL<-cluster.vcov(lm, data$id)
  coef<-coeftest(lm, vcovCL)
  acme_estimates<-coef %>% tidy() %>% filter(estimate!= "NA")%>% filter(term!="(Intercept)")%>% 
    mutate(lower=estimate-1.96*std.error) %>% mutate(upper=estimate+1.96*std.error) %>% 
    mutate(cate = sub(":.*", "", term)) %>% mutate(feature_level=NA)
  
  acme_estimates$feature_level[acme_estimates$term=="whotaxpaidYou would pay the tax to the Local Government Authority"]<-"Who\nLocal Government Authority"                                                     
  acme_estimates$feature_level[acme_estimates$term=="whotaxpaidYou would pay the tax to the State Internal Revenue Services"]<-"Who\nState Internal Revenue Services"                                              
  acme_estimates$feature_level[acme_estimates$term=="leveltaxThe tax on income would be between N24000 and N36,000"]<-"Level\nN24000-N36000"                                           
  acme_estimates$feature_level[acme_estimates$term=="leveltaxThe tax on income would be between N6000 and N12,000"]<-"Level\nN6000-N12000"                                                         
  acme_estimates$feature_level[acme_estimates$term=="deductYou would receive low-cost healthcare, in which you co-pay N 750 each time you go to the doctor"]<-"Deduct\nLow-cost healthcare"                
  acme_estimates$feature_level[acme_estimates$term=="howtaxYou would pay the tax to the bank"]<-"How\n bank"                                                                               
  acme_estimates$feature_level[acme_estimates$term=="howtaxYou would pay via UUSD, which uses telecom service just like when you load credit on your cell phone"]<-"How\n UUSD"            
  acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare already 3 months before you start paying income tax"]<-"Health Coverage\n3 months before you pay"                             
  #acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare at the moment you start paying income tax"]<-"Health Coverage\nat the moment start paying"                                       
  acme_estimates$feature_level[acme_estimates$term=="perkWhen paying income tax, you would receive a state identity card that you could access new government services with"]<-"Perk\nState identity card"
  acme_estimates$feature_level[acme_estimates$term=="perkYour health coverage would be extended to one partner and four children"]<-"Perk\nHealth coverage extended to 1 partner and 4 children"
  acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare 3 months after you start paying income tax"]<-"Health Coverage\n3 months after you pay"
  
  acme_estimates %>% 
    filter(estimate!=0) %>% ## different style: remove baseline
    ggplot(aes(x=feature_level, y=estimate, color=feature_level)) +
    #ggplot(aes(x=term, y=estimate)) +
    geom_point(position=position_dodge(width = .3), size=4) +
    geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .3), size=0.7) +
    geom_hline(yintercept=0, linetype="dashed", color="grey") +
    theme(legend.title=element_blank(),
          legend.text=element_text(size=15),
          axis.title=element_text(size=15),
          axis.text=element_text(size=15))+
    theme(axis.title.y = element_blank()) + theme(axis.title.x = element_blank()) +
    theme(legend.position = "none") +
    coord_flip() -> p
  
  return(acme_estimates)
}


benefitlist <-c("suffer_chronic_disease", "accessed_healthcare", "enrolled_insurance", "confidence_health_clinic_dummy","trust_health_institutions2")

for (j in benefitlist){
  name = deparse(substitute(j))
  data_var <- subset(data, get(j)==1)
  nam <-paste0("data_",j)
  assign(nam, data_var)
}    

for (j in benefitlist){
  name = deparse(substitute(j))
  data_var <- subset(data, get(j)==0)
  nam <-paste0("data_",j,"_no")
  assign(nam, data_var)
}    

for (i in benefitlist) {
  nam <-paste0("data_",i)
  dat <- get(nam)
  result<-acme_reg_old(dat)
  nam <-paste0(i,"_est")
  assign(nam, result)
}

for (i in benefitlist) {
  nam <-paste0("data_",i,"_no")
  dat <- get(nam)
  result<-acme_reg_old(dat)
  nam <-paste0(i,"_est_no")
  assign(nam, result)
}


suffer_chronic_disease_est_org <- suffer_chronic_disease_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
suffer_chronic_disease_est<-  suffer_chronic_disease_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
suffer_chronic_disease_est <- suffer_chronic_disease_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
suffer_chronic_disease_est <- suffer_chronic_disease_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
suffer_chronic_disease_est <- suffer_chronic_disease_est %>% filter(cate=="hea" | cate=="lev")
suffer_chronic_disease_est <- suffer_chronic_disease_est %>% mutate(var = "Yes") 
suffer_chronic_disease_est_no<-  suffer_chronic_disease_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
suffer_chronic_disease_est_no <- suffer_chronic_disease_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
suffer_chronic_disease_est_no <- suffer_chronic_disease_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
suffer_chronic_disease_est_no <- suffer_chronic_disease_est_no %>% filter(cate=="hea" | cate=="lev")
suffer_chronic_disease_est_no <- suffer_chronic_disease_est_no %>% mutate(var = "No") 
suffer_chronic_disease_all <-rbind(suffer_chronic_disease_est,suffer_chronic_disease_est_no)
suffer_chronic_disease_all <- suffer_chronic_disease_all %>% mutate(class = "suffer_chronic_disease_all")

accessed_healthcare_est_org<-  accessed_healthcare_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
accessed_healthcare_est<-  accessed_healthcare_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
accessed_healthcare_est <- accessed_healthcare_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
accessed_healthcare_est <- accessed_healthcare_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
accessed_healthcare_est <- accessed_healthcare_est %>% filter(cate=="hea" | cate=="lev")
accessed_healthcare_est <- accessed_healthcare_est %>% mutate(var = "Yes") 
accessed_healthcare_est_no<-  accessed_healthcare_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
accessed_healthcare_est_no <- accessed_healthcare_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
accessed_healthcare_est_no <- accessed_healthcare_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
accessed_healthcare_est_no <- accessed_healthcare_est_no %>% filter(cate=="hea" | cate=="lev")
accessed_healthcare_est_no <- accessed_healthcare_est_no %>% mutate(var = "No") 
accessed_healthcare_all <-rbind(accessed_healthcare_est,accessed_healthcare_est_no)
accessed_healthcare_all <- accessed_healthcare_all %>% mutate(class = "accessed_healthcare_all")

enrolled_insurance_est_org<-  enrolled_insurance_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
enrolled_insurance_est<-  enrolled_insurance_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
enrolled_insurance_est <- enrolled_insurance_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
enrolled_insurance_est <- enrolled_insurance_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
enrolled_insurance_est <- enrolled_insurance_est %>% filter(cate=="hea" | cate=="lev")
enrolled_insurance_est <- enrolled_insurance_est %>% mutate(var = "Yes") 
enrolled_insurance_est_no<-  enrolled_insurance_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
enrolled_insurance_est_no <- enrolled_insurance_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
enrolled_insurance_est_no <- enrolled_insurance_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
enrolled_insurance_est_no <- enrolled_insurance_est_no %>% filter(cate=="hea" | cate=="lev")
enrolled_insurance_est_no <- enrolled_insurance_est_no %>% mutate(var = "No") 
enrolled_insurance_all <-rbind(enrolled_insurance_est,enrolled_insurance_est_no)
enrolled_insurance_all <- enrolled_insurance_all %>% mutate(class = "enrolled_insurance_all")

confidence_health_clinic_dummy_est_org<-  confidence_health_clinic_dummy_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
confidence_health_clinic_dummy_est<-  confidence_health_clinic_dummy_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
confidence_health_clinic_dummy_est <- confidence_health_clinic_dummy_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
confidence_health_clinic_dummy_est <- confidence_health_clinic_dummy_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
confidence_health_clinic_dummy_est <- confidence_health_clinic_dummy_est %>% filter(cate=="hea" | cate=="lev")
confidence_health_clinic_dummy_est <- confidence_health_clinic_dummy_est %>% mutate(var = "Yes") 
confidence_health_clinic_dummy_est_no<-  confidence_health_clinic_dummy_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
confidence_health_clinic_dummy_est_no <- confidence_health_clinic_dummy_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
confidence_health_clinic_dummy_est_no <- confidence_health_clinic_dummy_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
confidence_health_clinic_dummy_est_no <- confidence_health_clinic_dummy_est_no %>% filter(cate=="hea" | cate=="lev")
confidence_health_clinic_dummy_est_no <- confidence_health_clinic_dummy_est_no %>% mutate(var = "No") 
confidence_health_clinic_dummy_all <-rbind(confidence_health_clinic_dummy_est,confidence_health_clinic_dummy_est_no)
confidence_health_clinic_dummy_all <- confidence_health_clinic_dummy_all %>% mutate(class = "confidence_health_clinic_dummy_all")

trust_health_institutions2_est_org<-  trust_health_institutions2_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
trust_health_institutions2_est<-  trust_health_institutions2_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
trust_health_institutions2_est <- trust_health_institutions2_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
trust_health_institutions2_est <- trust_health_institutions2_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
trust_health_institutions2_est <- trust_health_institutions2_est %>% filter(cate=="hea" | cate=="lev")
trust_health_institutions2_est <- trust_health_institutions2_est %>% mutate(var = "Yes") 
trust_health_institutions2_est_no<-  trust_health_institutions2_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
trust_health_institutions2_est_no <- trust_health_institutions2_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
trust_health_institutions2_est_no <- trust_health_institutions2_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
trust_health_institutions2_est_no <- trust_health_institutions2_est_no %>% filter(cate=="hea" | cate=="lev")
trust_health_institutions2_est_no <- trust_health_institutions2_est_no %>% mutate(var = "No") 
trust_health_institutions2_all <-rbind(trust_health_institutions2_est,trust_health_institutions2_est_no)
trust_health_institutions2_all <- trust_health_institutions2_all %>% mutate(class = "trust_health_institutions2_all")

benefit2<-rbind(suffer_chronic_disease_all, accessed_healthcare_all, enrolled_insurance_all, confidence_health_clinic_dummy_all,trust_health_institutions2_all)


benefit2 <- benefit2 %>% mutate(feature_level = factor(feature_level,levels=c( 
  "Health Coverage\n3 months before you pay","Health Coverage\n3 months after you pay","Health Coverage\nat the moment start paying",
  "Level\nN24000-N36000","Level\nN12000-N24000","Level\nN6000-N12000"
)))


benefit2_lev <- benefit2  %>% filter(cate=="lev")

class_label <-c(accessed_healthcare_all = "Access to Healthcare",suffer_chronic_disease_all = "Chronic Disease",enrolled_insurance_all = "Insurance", confidence_health_clinic_dummy_all = "Confidence in\nHealth Clinic",trust_health_institutions2_all = "Trust in\nHealth Institutions")

benefit2_lev %>% 
  ggplot(aes(estimate, feature_level, color=as.factor(var))) +
  #geom_point(position=position_dodge(width = .3),size=1) +
  ggstance::geom_pointrangeh(aes(xmin=lower, xmax=upper), position=position_dodge(width = .4), size=0.7) +
  geom_vline(xintercept=0, linetype="dashed", color="grey") +
  theme(legend.title=element_blank(),
        legend.text=element_text(size=20),
        axis.title=element_text(size=20),
        axis.text=element_text(size=20))+
  xlab("")+
  ylab("")+
  facet_grid(~ class, labeller = labeller (class = class_label))+
  theme(strip.text.x = element_text(size = 14, colour = "black"),
        axis.text.x = element_text(size = 14,face="bold"),
        axis.text.y = element_text(size = 20,face="bold")) + 
  guides(color=guide_legend(reverse = TRUE))-> p
p

####################################
############Table 3 ###############
####################################

acme_rounded <- benefit2_lev

numeric_cols <- c("estimate", "std.error", "statistic", "p.value", "lower", "upper")
for(col in numeric_cols) {
  acme_rounded[[col]] <- round(acme_rounded[[col]], 2)
}

print(acme_rounded)



####################################
############Figure 3 ###############
####################################



distancelist <-c("local_accountability", "service_satisfaction", "receive_social_benefits", "paid_tax", "political_engagement", "local_responsiveness")


for (j in distancelist){
  name = deparse(substitute(j))
  data_var <- subset(data, get(j)==1)
  nam <-paste0("data_",j)
  assign(nam, data_var)
}    

for (j in distancelist){
  name = deparse(substitute(j))
  data_var <- subset(data, get(j)==0)
  nam <-paste0("data_",j,"_no")
  assign(nam, data_var)
}    

for (i in distancelist) {
  nam <-paste0("data_",i)
  dat <- get(nam)
  result<-acme_reg_old(dat)
  nam <-paste0(i,"_est")
  assign(nam, result)
}

for (i in distancelist) {
  nam <-paste0("data_",i,"_no")
  dat <- get(nam)
  result<-acme_reg_old(dat)
  nam <-paste0(i,"_est_no")
  assign(nam, result)
}


service_satisfaction_est<-  service_satisfaction_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
service_satisfaction_est <- service_satisfaction_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
service_satisfaction_est <- service_satisfaction_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
service_satisfaction_est <- service_satisfaction_est %>% filter(cate=="hea" | cate=="lev")
service_satisfaction_est <- service_satisfaction_est %>% mutate(var = "Yes") 
service_satisfaction_est_no<-  service_satisfaction_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
service_satisfaction_est_no <- service_satisfaction_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
service_satisfaction_est_no <- service_satisfaction_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
service_satisfaction_est_no <- service_satisfaction_est_no %>% filter(cate=="hea" | cate=="lev")
service_satisfaction_est_no <- service_satisfaction_est_no %>% mutate(var = "No") 
service_satisfaction_all <-rbind(service_satisfaction_est,service_satisfaction_est_no)
service_satisfaction_all <- service_satisfaction_all %>% mutate(class = "service_satisfaction_all")

political_engagement_est<-  political_engagement_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
political_engagement_est <- political_engagement_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
political_engagement_est <- political_engagement_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
political_engagement_est <- political_engagement_est %>% filter(cate=="hea" | cate=="lev")
political_engagement_est <- political_engagement_est %>% mutate(var = "Yes") 
political_engagement_est_no<-  political_engagement_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
political_engagement_est_no <- political_engagement_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
political_engagement_est_no <- political_engagement_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
political_engagement_est_no <- political_engagement_est_no %>% filter(cate=="hea" | cate=="lev")
political_engagement_est_no <- political_engagement_est_no %>% mutate(var = "No") 
political_engagement_all <-rbind(political_engagement_est,political_engagement_est_no)
political_engagement_all <- political_engagement_all %>% mutate(class = "political_engagement_all")

receive_social_benefits_est<-  receive_social_benefits_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
receive_social_benefits_est <- receive_social_benefits_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
receive_social_benefits_est <- receive_social_benefits_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
receive_social_benefits_est <- receive_social_benefits_est %>% filter(cate=="hea" | cate=="lev")
receive_social_benefits_est <- receive_social_benefits_est %>% mutate(var = "Yes") 
receive_social_benefits_est_no<-  receive_social_benefits_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
receive_social_benefits_est_no <- receive_social_benefits_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
receive_social_benefits_est_no <- receive_social_benefits_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
receive_social_benefits_est_no <- receive_social_benefits_est_no %>% filter(cate=="hea" | cate=="lev")
receive_social_benefits_est_no <- receive_social_benefits_est_no %>% mutate(var = "No") 
receive_social_benefits_all <-rbind(receive_social_benefits_est,receive_social_benefits_est_no)
receive_social_benefits_all <- receive_social_benefits_all %>% mutate(class = "receive_social_benefits_all")


paid_tax_est<-  paid_tax_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
paid_tax_est <- paid_tax_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
paid_tax_est <- paid_tax_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
paid_tax_est <- paid_tax_est %>% filter(cate=="hea" | cate=="lev")
paid_tax_est <- paid_tax_est %>% mutate(var = "Yes") 
paid_tax_est_no<-  paid_tax_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
paid_tax_est_no <- paid_tax_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
paid_tax_est_no <- paid_tax_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
paid_tax_est_no <- paid_tax_est_no %>% filter(cate=="hea" | cate=="lev")
paid_tax_est_no <- paid_tax_est_no %>% mutate(var = "No") 
paid_tax_all <-rbind(paid_tax_est,paid_tax_est_no)
paid_tax_all <- paid_tax_all %>% mutate(class = "paid_tax_all")

local_responsiveness_est<-  local_responsiveness_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
local_responsiveness_est <- local_responsiveness_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
local_responsiveness_est <- local_responsiveness_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
local_responsiveness_est <- local_responsiveness_est %>% filter(cate=="hea" | cate=="lev")
local_responsiveness_est <- local_responsiveness_est %>% mutate(var = "Yes") 
local_responsiveness_est_no<-  local_responsiveness_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
local_responsiveness_est_no <- local_responsiveness_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
local_responsiveness_est_no <- local_responsiveness_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
local_responsiveness_est_no <- local_responsiveness_est_no %>% filter(cate=="hea" | cate=="lev")
local_responsiveness_est_no <- local_responsiveness_est_no %>% mutate(var = "No") 
local_responsiveness_all <-rbind(local_responsiveness_est,local_responsiveness_est_no)
local_responsiveness_all <- local_responsiveness_all %>% mutate(class = "local_responsiveness_all")

local_accountability_est<-  local_accountability_est %>% mutate(cate = substr(term, start = 1, stop = 3)) 
local_accountability_est <- local_accountability_est %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
local_accountability_est <- local_accountability_est %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
local_accountability_est <- local_accountability_est %>% filter(cate=="hea" | cate=="lev")
local_accountability_est <- local_accountability_est %>% mutate(var = "Yes") 
local_accountability_est_no<-  local_accountability_est_no %>% mutate(cate = substr(term, start = 1, stop = 3)) 
local_accountability_est_no <- local_accountability_est_no %>% add_row(feature_level = "Level\nN12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
local_accountability_est_no <- local_accountability_est_no %>% add_row(feature_level = "Health Coverage\nat the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
local_accountability_est_no <- local_accountability_est_no %>% filter(cate=="hea" | cate=="lev")
local_accountability_est_no <- local_accountability_est_no %>% mutate(var = "No") 
local_accountability_all <-rbind(local_accountability_est,local_accountability_est_no)
local_accountability_all <- local_accountability_all %>% mutate(class = "local_accountability_all")

distance<-rbind(local_accountability_all, service_satisfaction_all, receive_social_benefits_all, paid_tax_all, political_engagement_all, local_responsiveness_all)

distance <- distance %>% mutate(feature_level = factor(feature_level,levels=c( 
  "Health Coverage\n3 months before you pay","Health Coverage\n3 months after you pay","Health Coverage\nat the moment start paying",
  "Level\nN24000-N36000","Level\nN12000-N24000","Level\nN6000-N12000"
)))


distance_lev <- distance  %>% filter(cate=="lev")


class_label <-c(service_satisfaction_all = "Service Satisfaction", political_engagement_all = "Political engagement",local_responsiveness_all = "Local Responsiveness", receive_social_benefits_all = "Social Benefits",paid_tax_all = "Tax Payment", local_accountability_all = "Local Accountability")

distance_lev %>% 
  ggplot(aes(estimate, feature_level, color=as.factor(var))) +
  #geom_point(position=position_dodge(width = .3),size=1) +
  ggstance::geom_pointrangeh(aes(xmin=lower, xmax=upper), position=position_dodge(width = .4), size=0.7) +
  geom_vline(xintercept=0, linetype="dashed", color="grey") +
  theme(legend.title=element_blank(),
        legend.text=element_text(size=20),
        axis.title=element_text(size=20),
        axis.text=element_text(size=20))+
  xlab("")+
  ylab("")+
  facet_grid(~ class, labeller = labeller (class = class_label))+
  scale_x_continuous(breaks=c(-0.1, 0, 0.1))+
  theme(strip.text.x = element_text(size = 14, colour = "black"),
        axis.text.x = element_text(size = 14,face="bold"),
        axis.text.y = element_text(size = 20,face="bold")) + 
  guides(color=guide_legend(reverse = TRUE))-> p
p

####################################
############Table 4 ###############
####################################

acme_rounded <- distance_lev

numeric_cols <- c("estimate", "std.error", "statistic", "p.value", "lower", "upper")
for(col in numeric_cols) {
  acme_rounded[[col]] <- round(acme_rounded[[col]], 2)
}

print(acme_rounded)



####################################
############APPENDIX ###############
####################################

####################################
############Figure A10 ##############
####################################

lm<-lm(rating_bi ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data, weights=data$sampling_weight)
vcovCL<-cluster.vcov(lm, data$id)
coef<-coeftest(lm, vcovCL)

acme_estimates<-coef %>% tidy() %>% filter(term!="(Intercept)")%>% 
  mutate(lower=estimate-1.96*std.error) %>% mutate(upper=estimate+1.96*std.error) %>% 
  mutate(cate = substr(term, start = 1, stop = 3)) %>% mutate(feature_level=NA)

#  %>% filter(estimate!= "NA")
acme_estimates$feature_level[acme_estimates$term=="whotaxpaidYou would pay the tax to the Local Government Authority"]<-"Local Government Authority"                                                     
acme_estimates$feature_level[acme_estimates$term=="whotaxpaidYou would pay the tax to the State Internal Revenue Services"]<-"State Internal Revenue Services"                                              
acme_estimates$feature_level[acme_estimates$term=="leveltaxThe tax on income would be between N24000 and N36,000"]<-"N24000-N36000"                                           
acme_estimates$feature_level[acme_estimates$term=="leveltaxThe tax on income would be between N6000 and N12,000"]<-"N6000-N12000"                                                         
acme_estimates$feature_level[acme_estimates$term=="deductYou would receive low-cost healthcare, in which you co-pay N 750 each time you go to the doctor"]<-"Low-cost healthcare"                
acme_estimates$feature_level[acme_estimates$term=="howtaxYou would pay the tax to the bank"]<-"Bank"                                                                               
acme_estimates$feature_level[acme_estimates$term=="howtaxYou would pay via UUSD, which uses telecom service just like when you load credit on your cell phone"]<-"UUSD"            
acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare already 3 months before you start paying income tax"]<-"3 months before you pay"                             
#acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare at the moment you start paying income tax"]<-"at the moment start paying"                                       
acme_estimates$feature_level[acme_estimates$term=="perkWhen paying income tax, you would receive a state identity card that you could access new government services with"]<-"State identity card"
acme_estimates$feature_level[acme_estimates$term=="perkYour health coverage would be extended to one partner and four children"]<-"Health coverage extended to 1 partner and 4 children"
acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare 3 months after you start paying income tax"]<-"3 months after you pay"

acme_estimates <- acme_estimates %>% add_row(feature_level = "Who", estimate = NA, lower= 0, upper=0, cate="who")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Level", estimate = NA, lower= 0, upper=0, cate="lev")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Deduct", estimate = NA, lower= 0, upper=0, cate="ded")
acme_estimates <- acme_estimates %>% add_row(feature_level = "How", estimate = NA, lower= 0, upper=0, cate="how")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Health Coverage", estimate = NA, lower= 0, upper=0, cate="hea")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Perk", estimate = NA, lower= 0, upper=0, cate="per")

acme_estimates <- acme_estimates %>% add_row(feature_level = "Cooperative or association", estimate = 0, lower= 0, upper=0, cate="who")
acme_estimates <- acme_estimates %>% add_row(feature_level = "N12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Free healthcare", estimate = 0, lower= 0, upper=0, cate="ded")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Tax officer or collecting agent", estimate = 0, lower= 0, upper=0, cate="how")
acme_estimates <- acme_estimates %>% add_row(feature_level = "At the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
acme_estimates <- acme_estimates %>% add_row(feature_level = "Healthcare is kept for a certain time when income is lost", estimate = 0, lower= 0, upper=0, cate="per")


acme_estimates <- acme_estimates %>% mutate(feature_level = factor(feature_level,levels=c( 
  "State identity card", "Health coverage extended to 1 partner and 4 children", "Healthcare is kept for a certain time when income is lost","Perk",
  "3 months before you pay",  "3 months after you pay","At the moment start paying","Health Coverage",
  "Bank","UUSD", "Tax officer or collecting agent","How",
  "Low-cost healthcare","Free healthcare","Deduct",
  "N24000-N36000", "N6000-N12000" ,"N12000-N24000","Level",
  "Local Government Authority","State Internal Revenue Services","Cooperative or association","Who")))

acme_estimates %>% 
  #filter(estimate!=0) %>% 
  ggplot(aes(estimate, feature_level)) +
  geom_point(position=position_dodge(width = .4),size=4, alpha=0.8) +
  ggstance::geom_pointrangeh(aes(xmin=lower, xmax=upper), position=position_dodge(width = .4), size=0.7) +
  geom_vline(xintercept=0, linetype="dashed", color="grey") +
  theme(legend.title=element_blank(),
        legend.text=element_text(size=15),
        axis.title=element_text(size=15),
        axis.text=element_text(size=15))+
  xlab("")+
  ylab("")+
  scale_y_discrete(labels=rev(c(expression(paste(bold("Who:"))),
                                expression(paste("Cooperative or association")),
                                expression(paste("State Internal Revenue Services")),
                                expression(paste("Local Government Authority")),
                                expression(paste(bold("Level:"))),
                                expression(paste("N12000-N24000")),
                                expression(paste("N6000-N12000")),
                                expression(paste("N24000-N36000")),
                                expression(paste(bold("Deduct:"))),
                                expression(paste("Free healthcare")),
                                expression(paste("Low-cost healthcare")),
                                expression(paste(bold("How:"))),
                                expression(paste("Tax officer or collecting agent")),
                                expression(paste("UUSD")),
                                expression(paste("Bank")),
                                expression(paste(bold("Health Coverage:"))),
                                expression(paste("At the moment start paying")),
                                expression(paste("3 months after you pay")),
                                expression(paste( "3 months before you pay")),
                                expression(paste(bold("Perk:"))),
                                expression(paste("Healthcare is kept for a certain time when income is lost")),
                                expression(paste("Health coverage extended to 1 partner and 4 children")),
                                expression(paste("State identity card"))
  )))+
  guides(color=guide_legend(reverse = TRUE))-> p
p


####################################
############Figure A11 ##############
####################################

health2_lev <- benefit2  %>% filter(cate=="hea")

class_label <-c(accessed_healthcare_all = "Access to Healthcare",suffer_chronic_disease_all = "Chronic Disease",enrolled_insurance_all = "Insurance", confidence_health_clinic_dummy_all = "Confidence in\nHealth Clinic",trust_health_institutions2_all = "Trust in\nHealth Institutions")

health2_lev %>% 
  ggplot(aes(estimate, feature_level, color=as.factor(var))) +
  #geom_point(position=position_dodge(width = .3),size=1) +
  ggstance::geom_pointrangeh(aes(xmin=lower, xmax=upper), position=position_dodge(width = .4), size=0.7) +
  geom_vline(xintercept=0, linetype="dashed", color="grey") +
  theme(legend.title=element_blank(),
        legend.text=element_text(size=20),
        axis.title=element_text(size=20),
        axis.text=element_text(size=20))+
  xlab("")+
  ylab("")+
  facet_grid(~ class, labeller = labeller (class = class_label))+
  theme(strip.text.x = element_text(size = 14, colour = "black"),
        axis.text.x = element_text(size = 14,face="bold"),
        axis.text.y = element_text(size = 20,face="bold")) + 
  guides(color=guide_legend(reverse = TRUE))-> p
p


####################################
############Figure A12 ##############
####################################

distance_hea <- distance  %>% filter(cate=="hea")

class_label <-c(service_satisfaction_all = "Service Satisfaction", political_engagement_all = "Political engagement",local_responsiveness_all = "Local Responsiveness", receive_social_benefits_all = "Social Benefits",paid_tax_all = "Tax Payment", local_accountability_all = "Local Accountability")


distance_hea %>% 
  ggplot(aes(estimate, feature_level, color=as.factor(var))) +
  #geom_point(position=position_dodge(width = .3),size=1) +
  ggstance::geom_pointrangeh(aes(xmin=lower, xmax=upper), position=position_dodge(width = .4), size=0.7) +
  geom_vline(xintercept=0, linetype="dashed", color="grey") +
  theme(legend.title=element_blank(),
        legend.text=element_text(size=20),
        axis.title=element_text(size=20),
        axis.text=element_text(size=20))+
  xlab("")+
  ylab("")+
  scale_x_continuous(breaks=c(-0.1, 0, 0.1))+
  facet_grid(~ class, labeller = labeller (class = class_label))+
  theme(strip.text.x = element_text(size = 14, colour = "black"),
        axis.text.x = element_text(size = 14,face="bold"),
        axis.text.y = element_text(size = 20,face="bold")) + 
  guides(color=guide_legend(reverse = TRUE))-> p
p


####################################
############Figure A13-A15 ##############
####################################

acme_reg <- function(data_input){
  data<-data_input
  lm<-lm(rating ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data)
  vcovCL<-cluster.vcov(lm, data$id)
  coef<-coeftest(lm, vcovCL)
  acme_estimates<-coef %>% tidy() %>% filter(term!="(Intercept)")%>% 
    mutate(lower=estimate-1.96*std.error) %>% mutate(upper=estimate+1.96*std.error) %>% 
    mutate(cate = substr(term, start = 1, stop = 3)) %>% mutate(feature_level=NA)
  
  #  %>% filter(estimate!= "NA")
  acme_estimates$feature_level[acme_estimates$term=="whotaxpaidYou would pay the tax to the Local Government Authority"]<-"Local Government Authority"                                                     
  acme_estimates$feature_level[acme_estimates$term=="whotaxpaidYou would pay the tax to the State Internal Revenue Services"]<-"State Internal Revenue Services"                                              
  acme_estimates$feature_level[acme_estimates$term=="leveltaxThe tax on income would be between N24000 and N36,000"]<-"N24000-N36000"                                           
  acme_estimates$feature_level[acme_estimates$term=="leveltaxThe tax on income would be between N6000 and N12,000"]<-"N6000-N12000"                                                         
  acme_estimates$feature_level[acme_estimates$term=="deductYou would receive low-cost healthcare, in which you co-pay N 750 each time you go to the doctor"]<-"Low-cost healthcare"                
  acme_estimates$feature_level[acme_estimates$term=="howtaxYou would pay the tax to the bank"]<-"Bank"                                                                               
  acme_estimates$feature_level[acme_estimates$term=="howtaxYou would pay via UUSD, which uses telecom service just like when you load credit on your cell phone"]<-"UUSD"            
  acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare already 3 months before you start paying income tax"]<-"3 months before you pay"                             
  #acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare at the moment you start paying income tax"]<-"at the moment start paying"                                       
  acme_estimates$feature_level[acme_estimates$term=="perkWhen paying income tax, you would receive a state identity card that you could access new government services with"]<-"State identity card"
  acme_estimates$feature_level[acme_estimates$term=="perkYour health coverage would be extended to one partner and four children"]<-"Health coverage extended to 1 partner and 4 children"
  acme_estimates$feature_level[acme_estimates$term=="healthcovYou would receive healthcare 3 months after you start paying income tax"]<-"3 months after you pay"
  
  acme_estimates <- acme_estimates %>% add_row(feature_level = "Who", estimate = NA, lower= 0, upper=0, cate="who")
  acme_estimates <- acme_estimates %>% add_row(feature_level = "Level", estimate = NA, lower= 0, upper=0, cate="lev")
  acme_estimates <- acme_estimates %>% add_row(feature_level = "Deduct", estimate = NA, lower= 0, upper=0, cate="ded")
  acme_estimates <- acme_estimates %>% add_row(feature_level = "How", estimate = NA, lower= 0, upper=0, cate="how")
  acme_estimates <- acme_estimates %>% add_row(feature_level = "Health Coverage", estimate = NA, lower= 0, upper=0, cate="hea")
  acme_estimates <- acme_estimates %>% add_row(feature_level = "Perk", estimate = NA, lower= 0, upper=0, cate="per")
  
  acme_estimates <- acme_estimates %>% add_row(feature_level = "Cooperative or association", estimate = 0, lower= 0, upper=0, cate="who")
  acme_estimates <- acme_estimates %>% add_row(feature_level = "N12000-N24000", estimate = 0, lower= 0, upper=0, cate="lev")
  acme_estimates <- acme_estimates %>% add_row(feature_level = "Free healthcare", estimate = 0, lower= 0, upper=0, cate="ded")
  acme_estimates <- acme_estimates %>% add_row(feature_level = "Tax officer or collecting agent", estimate = 0, lower= 0, upper=0, cate="how")
  acme_estimates <- acme_estimates %>% add_row(feature_level = "At the moment start paying", estimate = 0, lower= 0, upper=0, cate="hea")
  acme_estimates <- acme_estimates %>% add_row(feature_level = "Healthcare is kept for a certain time when income is lost", estimate = 0, lower= 0, upper=0, cate="per")
  
  
  acme_estimates <- acme_estimates %>% mutate(feature_level = factor(feature_level,levels=c( 
    "State identity card", "Health coverage extended to 1 partner and 4 children", "Healthcare is kept for a certain time when income is lost","Perk",
    "3 months before you pay",  "3 months after you pay","At the moment start paying","Health Coverage",
    "Bank","UUSD", "Tax officer or collecting agent","How",
    "Low-cost healthcare","Free healthcare","Deduct",
    "N24000-N36000", "N6000-N12000" ,"N12000-N24000","Level",
    "Local Government Authority","State Internal Revenue Services","Cooperative or association","Who")))
  
  return(acme_estimates)
}

draw_acme_hetero <- function(data_input,fname){
  all <- data_input
  all %>% 
    #filter(estimate!=0) %>% 
    ggplot(aes(estimate, feature_level, color=as.factor(var))) +
    geom_point(position=position_dodge(width = .4),size=4, alpha=0.8) +
    ggstance::geom_pointrangeh(aes(xmin=lower, xmax=upper), position=position_dodge(width = .4), size=0.7) +
    geom_vline(xintercept=0, linetype="dashed", color="grey") +
    theme(legend.title=element_blank(),
          legend.text=element_text(size=15),
          axis.title=element_text(size=15),
          axis.text=element_text(size=15))+
    xlab("")+
    ylab("")+
    scale_y_discrete(labels=rev(c(expression(paste(bold("Who:"))),
                                  expression(paste("Cooperative or association")),
                                  expression(paste("State Internal Revenue Services")),
                                  expression(paste("Local Government Authority")),
                                  expression(paste(bold("Level:"))),
                                  expression(paste("N12000-N24000")),
                                  expression(paste("N6000-N12000")),
                                  expression(paste("N24000-N36000")),
                                  expression(paste(bold("Deduct:"))),
                                  expression(paste("Free healthcare")),
                                  expression(paste("Low-cost healthcare")),
                                  expression(paste(bold("How:"))),
                                  expression(paste("Tax officer or collecting agent")),
                                  expression(paste("UUSD")),
                                  expression(paste("Bank")),
                                  expression(paste(bold("Health Coverage:"))),
                                  expression(paste("At the moment start paying")),
                                  expression(paste("3 months after you pay")),
                                  expression(paste( "3 months before you pay")),
                                  expression(paste(bold("Perk:"))),
                                  expression(paste("Healthcare is kept for a certain time when income is lost")),
                                  expression(paste("Health coverage extended to 1 partner and 4 children")),
                                  expression(paste("State identity card"))
    )))+
    guides(color=guide_legend(reverse = TRUE))-> p
  ggsave(p, filename=paste0("acme_hetero", fname, ".pdf"), 
         family="sans", width=30, height=30, units= "cm", dpi=800)
}

data_urban<-data[which(data$urban==1),]
data_rural<-data[which(data$urban==0),]

data_female<-data[which(data$gender=="Female"),]
data_male<-data[which(data$gender=="Male"),]

data_high_educ <-data[which(data$high_educ==1),]
data_mid_educ <-data[which(data$high_educ!=1 & data$low_educ!=1),]
data_low_educ <-data[which(data$low_educ==1),]


data_urban_est<-acme_reg(data_urban)
data_rural_est<-acme_reg(data_rural)

data_urban_est <- data_urban_est %>% mutate(var = "Urban") 
data_rural_est <- data_rural_est %>% mutate(var = "Rural") 

data_female_est<-acme_reg(data_female)
data_male_est<-acme_reg(data_male)

data_female_est <- data_female_est %>% mutate(var = "Female") 
data_male_est <- data_male_est %>% mutate(var = "Male") 

data_high_educ_est<-acme_reg(data_high_educ)
data_mid_educ_est<-acme_reg(data_mid_educ)
data_low_educ_est<-acme_reg(data_low_educ)

data_high_educ_est <- data_high_educ_est %>% mutate(var = "High education") 
data_mid_educ_est <- data_mid_educ_est %>% mutate(var = "Mid education") 
data_low_educ_est <- data_low_educ_est %>% mutate(var = "Low education") 


gender_all <-rbind(data_female_est,data_male_est)
draw_acme_hetero(gender_all,"gender")

urban_all <-rbind(data_urban_est,data_rural_est)
draw_acme_hetero(urban_all,"urban")

educ_all <-rbind(data_high_educ_est,data_mid_educ_est,data_low_educ_est)
draw_acme_hetero(educ_all,"educ")




####################################
############Table A2 ##############
####################################

data$perk <- factor(data$perk)
data$head_of_household <- factor(data$head_of_household)
data$gender <- factor(data$gender)
data$urban <- factor(data$urban)
data$schooling <- factor(data$schooling)
data$enrolled_insurance <- factor(data$enrolled_insurance)
data$medical_need <- factor(data$medical_need)
data$religion <- factor(data$religion)



lm<-lm(as.numeric(head_of_household) ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data, weights=data$sampling_weight)
vcovCL<-cluster.vcov(lm, data$id)
coef2<-coeftest(lm, vcovCL)

lm<-lm(as.numeric(gender) ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data, weights=data$sampling_weight)
vcovCL<-cluster.vcov(lm, data$id)
coef3<-coeftest(lm, vcovCL)

lm<-lm(as.numeric(urban) ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data, weights=data$sampling_weight)
vcovCL<-cluster.vcov(lm, data$id)
coef4<-coeftest(lm, vcovCL)

lm<-lm(as.numeric(schooling) ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data, weights=data$sampling_weight)
vcovCL<-cluster.vcov(lm, data$id)
coef5<-coeftest(lm, vcovCL)

lm<-lm(as.numeric(enrolled_insurance) ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data, weights=data$sampling_weight)
vcovCL<-cluster.vcov(lm, data$id)
coef6<-coeftest(lm, vcovCL)

lm<-lm(as.numeric(medical_need) ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data, weights=data$sampling_weight)
vcovCL<-cluster.vcov(lm, data$id)
coef7<-coeftest(lm, vcovCL)

lm<-lm(as.numeric(religion) ~ whotaxpaid + leveltax + deduct + howtax + healthcov + perk, data=data, weights=data$sampling_weight)
vcovCL<-cluster.vcov(lm, data$id)
coef8<-coeftest(lm, vcovCL)


stargazer(coef2,coef3,coef4,coef5,coef6,coef7,coef8,
          column.labels = c("Head of household","Female", "Urban", "Education","Enrolled insurance ","Medical need", "Religion"),
          float=F,
          keep.stat = "n",
          star.cutoffs = c(0.05, 0.01),
          column.sep.width = "5pt",
          notes.align = "l",
          out = "reg_balance.tex")
